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We present a numerical scheme for efficiently extracting the higher-order moments and cumulants 
of various operators on spin systems represented as tensor product states, for both finite and infinite 
systems, and present several applications for such quantities. For example, the second cumulant of 
the energy of a state, {AH^), gives a straightforward method to check the convergence of numerical 
ground-state approximation algorithms. Additionally, we discuss the use of moments and cumulants 
in the study of phase transitions. Of particular interest is the application of our method to calculate 
the so-called Binders cumulant, which we use to detect critical points and study the critical exponent 
of the correlation length with only small finite numerical calculations. We apply these methods to 
study the behavior of a family of one-dimensional models (the transverse Ising model, the spin-1 
Ising model, and the spin-1 Ising model in a crystal field), as well as the two-dimensional Ising 
model on a square lattice. Our results show that in one dimension, cumulant-based methods can 
produce precise estimates of the critical points at a low computational cost, and show promise for 
two-dimensional systems as well. 


I. INTRODUCTION 

The understanding of quantum many-body systems is 
one of the foremost goals of modern quantum physics. In 
addition to various analytical approaches, a new set of 
numerical tools has emerged for this purpose in recent 
years, based on tensor network representations of such 
systems m- These techniques take advantage of the 
so-called “area law” for entanglement entropy, obeyed 
by the low-energy eigenstates of gapped Hamiltonians 
with local interactions. Tensor networks naturally em¬ 
body this entanglement structure, dramatically simpli¬ 
fying the degrees of freedom required to describe such 
states. Initially introduced in the context of gapped one¬ 
dimensional systems, where they are referred to as “ma¬ 
trix product states” (MPS) ISHlj, tensor network meth¬ 
ods rose to even greater prominence when it was real¬ 
ized that the celebrated “Density Matrix Renormaliza¬ 
tion Group” (DMRG) technique [H [9] could be reformu¬ 
lated in terms of an MPS [lOHTSj . Additional algorithms 
soon followed DMRG, such as ‘Time Evolving Block Dec¬ 
imation” (TEBD) [H], as well as the infinite-system ana¬ 
logues “iDMRG” and “iTEBD” [IH [K]. Since then, 
tensor network methods have also been readily applied 
to critical systems be means of the “Multiscale Entan¬ 
glement Renormalization Anzatz” (MERA) [T71[TH], as 
well as to higher-dimensional systems, where they are 
generally termed “tensor network states,” or “projected- 
entangled pair states” (PEPS) [TM25] . A principle goal 
of these algorithms is to obtain precise approximations 
to ground-state wave functions, which can be used for 
many purposes. For example, one can compute various 
quantities and observables in an effort to detect phase 
transitions, a central problem in many-body physics Ei- 

Eo]. 

In the Landau symmetry-breaking paradigm for phase 
transitions, one first looks for an operator M whose ex¬ 
pectation value {M) can serve as the order parameter. 


i.e. a quantity whose behavior changes sharply across a 
critical point. When this order parameter is represented 
by a local operator, it can be computed efficiently on a 
tensor network state [2]. But while an expectation value 
is the most straightforward piece of information associ¬ 
ated with an operator and a state, there is considerably 
more information available which one may want to com¬ 
pute. For instance, one may wish to study the higher 
moments of the operator, /i„ = (M"). A related set of 
quantities called “cumulants,” typically labelled is 
also frequently of interest. An obvious example is the 
variance of the operator (AM^), which is simply the sec¬ 
ond cumulant K 2 = p 2 — Ti- Even more important to 
the search for phase transitions is the so-called “Binder 
cumulant,” first introduced by Kurt Binder in 1981 in a 
study of the classical Ising Model [31]. In many settings, 
such as thermal or disordered systems, it is considered to 
be one of the most accurate and reliable means of detect¬ 
ing a critical point , and it has since been applied 

to a wide variety of models [35VI42] . 

Computing these higher order moments and cumu¬ 
lants, however, is less straightforward. Direct calculation 
quickly becomes impractical for large n, since the num¬ 
ber of terms to evaluate can be exponential in n. In a 
classical system with a Hamiltonian Hq, one might define 
H{X) = Hq + XM, and relate the higher moments of M to 
the derivatives of an associated partition function, using 
(M") = (/3j^)"Tr(e“^'^*^'’‘^). In quantum systems, how¬ 
ever, this equation only holds when \Hq,M\ = 0, which 
is not true for a wide variety of physically interesting 
cases. Because of these barriers to direct calculation, us¬ 
age of techniques such as the powerful Binder cumulant 
has in the past been generally confined to studies based 
on quantum Monte Carlo [55] . 

The question naturally arises whether these quantities 
can be efficiently and systematically evaluated using the 
elegant structure of a tensor network state. Here, we 
demonstrate that the answer is yes. The feasibility of us- 
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ing matrix product states for computing the second mo¬ 
ment of Hamiltonians was already pointed out by 1. Mc¬ 
Culloch in the context of DMRG [15] , via the technique 
of so-called “matrix product operators” (MPO) [50]. In 
this work, we propose a simple and efficient method to al¬ 
low all general moments and cumulants to be evaluated 
for tensor network states, based on moment-generating 
and cumulant-generating functions. We demonstrate the 
calculation of moments and cumulants for finite one¬ 
dimensional states, and show that the method can also 
be used for per-site cumulants in the case of an infi¬ 
nite system. We also show how the techniques naturally 
generalize to finite systems in higher dimensions. These 
methods have a variety of useful applications which are 
demonstrated at length, including the use of the Binder 
and other cumulants to detect critical points to relatively 
high precision at a low numerical cost. We also apply 
the second cumulant of the energy to examine the con¬ 
vergence of numerical methods based on imaginary time 
evolution. 

This paper is organized as follows: In Section II we 
review moments and cumulants, in particular presenting 
the Binder cumulant and some of its applications. In 
Section III, we very briefly review the MBS formalism 
and discuss how to compute certain simple expectation 
values. Section IV demonstrates how to use these ex¬ 
pectation values to efficiently compute the moments and 
cumulants of general operators on an MBS. Section V 
contains examples of the method as applied to three dif¬ 
ferent spin-chain models (the transverse Ising model, the 
spin-1 Ising model, and the spin-1 Ising model in a crystal 
field), as well as a demonstration of the method as ap¬ 
plied to a two dimensional system (the transverse Ising 
model on a square lattice). Our results are summarized 
in Section VI. 


II. MOMENTS, CUMULANTS, AND THE 
BINDER CUMULANT 


A state lip) and an operator M collectively imply a 
probability distribution: the probability density function 
of ip in M-space. The expectation value (M) specifies the 
central value of the distribution, while the complete set 
of “Moments” defines the entire the shape [Hj. The 
moment of the distribution is defined to be /i„ = (M"); 
the first moment fj,i is the expectation value (M) itself. 

The cumulants of the distribution, form an alter¬ 
native but equivalent way of specifying its shape. These 
cumulants contain, in total, the same information as the 
moments; a complete set of either moments or cumulants 
completely specifies the distribution. Indeed, the cu¬ 
mulant can always be expressed as a polynomial combi¬ 
nation of the first n moments, and vice versa |45j . For 
example, as we have noted above, the second cumulant 
of the distribution, is the distribution’s variance, defined 
by 


K2 = M2 - Ml- (1) 

The third cumulant K 3 gives the distribution’s skew¬ 
ness, and is related to the first three moments by 

^^3 = M3 ~ 3m2Mi + 2mi- (2) 

Similarly, the fourth cumulant K 4 is related to the kur- 
tosis, and is given by 


Ki = 4/13/11 - 3m2 + 12m2Mi - 6Mi- (3) 

Although moments and cumulants are properly defined 
with respect to a distribution and hence depend on both 
M and \tp}, when \tp) is general or clear from context we 
shall refer to Mn («„) as the moment (cumulant) of 
M”. 


A. Binder’s Cumulant 

The aforementioned Binder Cumulant is a particularly 
useful quantity in the study of critical points and phase 
transitions. For some system with some known order 
parameter M, for example a total magnetization aj 
or a staggered magnetization Binder’s cu¬ 

mulant represents a modified version of that parameter’s 
4th cumulant. Though some slight variations exist in the 
definition, generally, it is given by 


C/4 = 1 - 


(M*) 

3(M2)2- 


(4) 


The utility of the Binder cumulant arises from the 
special features of its length dependance. The behav¬ 
ior of the Binder cumulant at a critical point depends 
only weakly on the size of the system, and elsewhere, 
its behavior with respect to the system size differs de¬ 
pending upon the phase. For example, below the critical 
point in a symmetry-breaking magnetic phase the cu¬ 
mulant will increase with the length of the system, but 
above the critical point, with symmetry unbroken, it de¬ 
creases instead. The result is that, when curves of the 
Binder Cumulant vs temperature are plotted for various 
lengths, the critical point is indicated by a simultaneous 
crossing. Typically, because the behavior at the critical 
point is already approximately universal, only a set of 
relatively small system sizes need be considered, elimi¬ 
nating the need for complicated extrapolations of very 
large systems to the thermodynamic limit. 

The Binder cumulant also gives access to the critical 
exponent of the correlation length, by means of tradi¬ 
tional finite size scaling techniques in which one seeks to 
“collapse” the data. mM- Up to some small finite 
size corrections (which become increasingly suppressed 
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as the system size increases), the cumulants show a stan¬ 
dard functional form, [46] 


Ui{L,B) = u(L^I''{B-BS), (5) 


3 

the expectation value (Q) is given by 


(7) 


where v is the usual critical exponent. A plot of 174 vs 
L^/'^{B—Bc) should therefore appear essentially indepen¬ 
dent of L, since all of the length-dependance has been ab¬ 
sorbed into the independent variable of the plot. Hence, 
Be and v can be treated as free parameters, and varied 
until this length-independence is optimized; for example, 
one could seek to minimize the total absolute square dis¬ 
tance between the curves for a variety of lengths L. 


III. MATRIX PRODUCT STATES AND 
EXPECTATION VALUES 


A. Finite-length chains 


(V’lQIV') = 


Tr 




or equivalently. 


( 8 ) 


= Tr 


L 

HE 

i=i sj.s' 


A 








(9) 


From this expression, it is clear that, up to normaliza¬ 
tion, the expectation value is simply a trace over a set of 
L “transfer matrices”; i.e. 


To demonstrate how to efficiently compute quantities 
such as the Binder cumulant for a system represented by 
an MPS, we must first review the nature of the MPS rep¬ 
resentation itself. Consider first a simple 1-dimensional 
spin chain of length L with periodic boundary conditions. 
As an MPS, this state will be expressed as 


(Q) 





where the T’s are defined as 


( 10 ) 


IV^) = (6) 

S 



( 11 ) 


In other words, a rank-three tensor “A^” has been as¬ 
sociated with each site “j”. One index of this tensor, de¬ 
noted Sj above, remains free and represents the physical 
degrees of freedom at site j. The other two indices, typ¬ 
ically termed the “virtual indices”, are contracted with 
the virtual indices of the neighboring tensors Aj_i and 
Aj_|_i. The dimension of these virtual indices, often la¬ 
belled X, is called the “bond dimension” of the MPS. 
The choice of x represents a numerical parameter which 
can be adjusted to suit the requirements of the context: 
smaller values are of course less numerically expensive, 
but larger values can allow the MPS to more accurately 
represent the features of the state, particularly in systems 
with long correlation lengths. 

In the one-dimensional case, since A is rank-3, for 
any fixed value of Sj the two remaining virtual indices 
simply represent a matrix. For this reason, the labels 
of the virtual indices are often suppressed, as in Eq. 
with the contraction represented by the expression 
in the same manner that one would 
write an ordinary product of matrices. 

Let us now consider how to calculate the expectation 
value of an operator with respect to a matrix product 
state. Consider first the case of a simple operator which 
is given by a tensor product of on-site operations. For 
such an operator, of the form 


This procedure is also demonstrated in graphical no¬ 
tation in Fig. [T] The norm of the state can be fixed in 
a similar fashion, by evaluating a transfer matrix for the 
special case where Qj =1. 

Tensor products of few-body operators can be handled 
in a similar fashion by grouping the relevant sites. More 
general operators are simply evaluated by decomposing 
them into a sum of tensor products. Considerably more 
detail on the general process of taking expectation values 
can be found in the now-extensive body of literature on 
matrix product states, [111117]. For our purposes, how¬ 
ever, it will be sufficient to be able to evaluate operators 
of the simple form in Eq. Q. 

B. Infinite-length chains 

One significant advantage of tensor network algorithms 
is how easily they allow one to directly study certain in¬ 
finite systems. This can be done for systems with some 
form of translation invariance, which can be completely 
represented by their unit cells. Eor example, consider 
an one-dimensional infinite system possessing translation 
invariance with respect to a unit cell of length i. Rep¬ 
resented as a matrix product, the state is of the form 
given by Eq. Q, but with the further restriction that 
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with a unit cell of length £, we shall consider only Q with 
Qj — Qj+i- 

In this situation, one can still sensibly define the ex¬ 
pectation value as 


(Q) 


(V'lV’) 



(14) 


Here, the transfer matrix is now “enlarged” to rep¬ 
resent an entire unit cell of the chain 


T, = \{T,. (15) 


FIG. 1. Graphical notation demonstrating the structures of 
matrix product states. In this notation, a shape represents a 
tensor, and a line represents an index. Connected lines be¬ 
tween shapes represent contracted indices between tensors. 

(a) A finite spin chain state li/i) represented as a matrix prod¬ 
uct state. The state is specified by the set of rank-three ten¬ 
sors {Aj}, with the physical degrees of freedom Sj left open. 

(b) The expectation value of a product operator Q = ®jQj 
with respect to |i/>). Each Qj acts locally on only one site. 
The total expectation value can be thought of as a trace over 
a product of transfer matrices Tj, defined in Eq. ( |11[ ). An 
example of an individual transfer matrix, T\ is highlighted. 


not all tensors Ai are distinct, and instead repeat ev¬ 
ery i sites. A state with only a one-site unit cell (full 
translation invariance) can therefore be specified by only 
a single tensor A 


In order to approach the infinite case, we shall first 
examine the case of a finite but very long chain of length 
L, so that the product in Eq. (14) is limited to L|^ terms, 

i.e. 


(Q)l = 



(16) 


For L sufficiently large, the product can then be ap¬ 
proximated by considering an eigenvalue decomposition 
of the transfer matrix Ti = UAU~^ and inserting it into 
Eq. ( |16[ ) (note that, by construction, the transfer matrix 
as defined by Eqs. (11) and (15) is Hermitian and hence 


diagonalizable). By applying the cyclic property of the 
trace operation, all U matrices can be made to cancel, 
leaving us only 


\tp) = |siS2...). 


( 12 ) 


Similarly, a state with two-site translation invariance 
' = 2) specified by two tensors, Ai, A 2 and has the form 


1^) =^Tr(A^)A(*=)A(^®)A^)...)|siS2S3S4...). (13) 


Of course, because the sums in Eqs. (12) and (13) run 


over an infinite number of sites, in general there is no way 
to specify or compute the coefficients. Certain expecta¬ 
tion values, on the other hand, may still be expressed as 
the limit of an infinite product of transfer matrices. It is 
quite common, for example, to consider the expectation 
value of an operator with the same translational invari¬ 
ance as the state in question, by looking at the per-site 
behavior. For our purposes, we will again be concerned 
with product operators of the form given in equation 
However, we will now restrict ourselves further by impos¬ 
ing translation invariance on Q. For an infinite system 




(aV<) , 


( 17 ) 


or, 


{Q)l = 


(V’lV’) 


Ed 

i=i 


L/l 


(18) 


where Xj are the diagonal elements of A, i.e. the eigenval¬ 
ues of the matrix . At this point, it can be observed that 
in the infinite limit, only the largest eigenvalue Xmax will 
contribute to the sum. In other words, we have simply 




(19) 


To fix the norm, we consider a pa rticular transfer ma¬ 
trix T), defined as usual by Eq. (11) for the special case 
where Q = 0^ Ij. Then, we calculate Xmax^ tbe largest 
eigenvalue of Tf, which satisfies 
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(V'lV') = (^Amaa;j ■ (20) 

Substituting, we have 

/A 

(Q)l=(^) . (21) 

We then gain access to the per-site behavior by means of 
a logarithm, which gives 


7 log(Q)L = log ( ^ 
^ V Aii 


l/£\ 


( 22 ) 


Eq. (231, however, does not depend on having a finite 


L. Thus, even for our infinite system, we can consider 
the limit 


lim j log{Q) = log (^ . (23) 

^ \ Xmax J 

Hence, with these procedures (illustrated graphically 
in Fig. [^, we can extract information about product op¬ 
erators in their infinite limit even though their expec¬ 
tation values generally diverge. As we will show below, 
this information will be sufficient to compute the cumu- 
lants of operators with translation invariance even in the 
infinite case. 


IV. EVALUATING HIGHER-ORDER 
MOMENTS AND CUMULANTS 



T£\R) = KM T<\R) = kr...\R') 


FIG. 2. (a) An infinite spin chain state possessing 

translation invariance with respect to a unit cell of length 
^ = 2, represented as a matrix product state, (b) A prod¬ 
uct operator Q = ®jQj which possesses the same translation 
symmetry as |'!/>); i.e. Qj = Qj+£ (c) To compute the quantity 
of interest, we first construct Ti, the transfer matrix contain¬ 
ing an entire unit cell of |'!/>) and Q, and extract its dominant 
eigenvalue Xmax- (d) To normalize the result, we will also 
need Tt (a transfer matrix which contains only the identity 
operator) and it’s dominant eigenvalue Xmax- The desired 
quantity limL_>oo ^ log(Q) is given by log(Xmax/Xmaxy^‘- 


From the result, it is clear that every (non-vanishing) 
moment (M") will appear in the expansion. Further¬ 
more, these moments can be directly accessed by com¬ 
puting 


A. Moment and Cumulant-Generating Functions 


F^^\a)=^^n+0{a), (25) 


For a given operator M and a state ji/i), there is an 
associated function F which contains all of the non-local 
information about the higher moments (M”), and yet, as 
we shall subsequently demonstrate, can still be efficiently 
evaluated within the framework of a tensor product state. 
In particular, this function is given by F{a) = (e“^). 

In probability theory, F{a) is termed the “moment gen¬ 
erating function” of the probability distribution. It is so 
named because the information about every moment of 
the distribution is not only contained, but readily acces¬ 
sible from this single function. This can be made explicit 
by considering a Taylor-expansion of about a = 0 
and then computing the expectation value in F{a) term- 
by-term 


F{a) 


1 -b a{M) 



(24) 


where (a) is as usual the derivative of F. 

The moment-generating function F is closely related to 
the so-called “characteristic function” of the distribution, 
G(a) = (e“^) [48]. For typical states with well-behaved 
wave functions, these functions will be essentially inter¬ 
changeable (up to a factor of i). Hence in this work, 
both will be used, sometimes in combination, depending 
on the particular moment or operator being computed. 
It should be noted, however, that for some “pathologi¬ 
cal” wave functions, such as those specifying a Lorentzian 
probability distribution, the function F(a) may fail to ex¬ 
ist. The characteristic function G(a), however, being the 
expectation value of a bounded operator, does not suffer 
from this complication in any situation |49] . 

While Eq. (25) gives a result for the desired moment 
which is only accurate up to first order in the parameter 
a, the precision can be improved by instead computing 
appropriate combinations of the functions F(a), F(—a). 
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G'(a), and G{—a). For example, when seeking to com¬ 
pute an even-ordered moment; i.e. a moment of the form 
we construct 


F{a)=F{a)+Fi-a) = 2 + a^{M^) + ^{M^) + ... (26) 

And hence obtain the desired moments from the rela¬ 
tion 


(X (M^")-tC>(a2). (27) 

Odd-ordered moments can of course be found to higher 
precision from F(a) — F(—a). Even greater precision can 
also be obtained by including the characteristic functions. 
The combination F{a) + F{—a) — G{a) — G{—a), for in¬ 
stance, determines (M^) up to 0(a'^). 

We employ a similar technique to extract the cumu- 
lants of the distribution. This is done by means of the 
“cumulant generating function,” defined as 


lp{a) = log F(a). (28) 

To see the utility of this function, observe that 

hia) « log(l + a{M) + ^(M^) + ...). (29) 

For small enough values of a, one can see from the 
expansion log(l + x) x — -I- ... that 

hia) = a{M) + ^(M)2 + ^(M2) + 0{a^). (30) 

Grouping these terms by the powers of a, we find that 
in fact 


Z_F(a) = OKI -I- —K2 + ... 


(31) 


In other words, the derivatives of the function Ip (a) 
give us direct access to the cumulants in the same manner 
as the moments in Eq. (251 


lp\a) = K„ + 0(a). (32) 

Of course, as with the moments, appropriate combina¬ 
tions of Ipict), lp{—a) and the associated complex func¬ 
tions can be used to suppress the higher order terms and 
improve the accuracy. 


B. Evaluating generating functions on a finite 
matrix product state 

Evaluating all of these moments and cumulants thus 
boils down to evaluating the expectation values of oper¬ 
ators like e“^. It remains to be shown that these opera¬ 
tors, which we term “moment generating operators,” can 
be applied in an efficient manner. Fortunately, the expo¬ 
nential structure of the operator guarantees that this is 
indeed the case. 

We will consider these operators in two cases, depend¬ 
ing on the nature of the operator M. The first, special 
case is the large class of operators where M can be writ¬ 
ten as Oj , where j runs over all the sites in the sys¬ 
tem and for some arbitrary set of on-site operators {Oj}. 
Most usefully, this set of operators contains the tradi¬ 
tional magnetization operators such as 
well as staggered magnetizations, crystal field magneti¬ 
zations, etc. Subsequently, we will examine the more 
general case where the terms within M act on more than 
one site at a time. 


1. Sums of single-body operators 

In this case, since the operators Oj all act at separate 
sites, the combined operator M is in fact simply a Kro- 
necker sum, M = 0^ Oj . From this it follows that we 
can write m 


aM 




(33) 


In other words, the moment-generating operator can 
be decomposed into a set of operators each acting 

only at a single site. The moment generating function, in 
turn, is just the expectation value of all these operators 
applied simultaneously. 

Since our moment-generating operator can be written 
in this tensor product format, we can compute its expec¬ 
tation value directly through Eqs. (10) and (II), with 
Qj = (see Fig. § . The calculation of our moment¬ 
generating function F(a) has therefore been reduced to 
the calculation of L transfer matrices and a single trace 
over their product. In practice, to calculate a higher 
moment like (M^), we then need to repeat this proce¬ 
dure and compute F for slightly different values of a, so 
that it is possible to evaluate the necessary derivative 
numerically. This can be done through any of the wide 
variety of standard methods; in this work we have used 
primarily the classic divided difference formulas [51) . In 
general, the more values of a at which F{a) is computed, 
the higher the accuracy of the derivative. However, since 
the initial a is already chosen to be quite small, in prac¬ 
tice it is often the case that only a very small number of 
points needs to be computed (to check the behavior and 
accuracy, the procedure can always be repeated with a 
smaller value of a). 
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Let us examine now the performance of the method 
for the case of (M^). A second derivative is necessary, 
which can be computed to second order in a from three 
values of a, centered at a = 0 m- Noting that when 
a = 0 we have trivially F{a) = 1, it follows that we only 
have to compute two expectation values, each involving 
the construction of (at most) L transfer matrices, which 
are then multiplied together and traced over. By com¬ 
parison, to compute (M^) directly, as the sum of all cor¬ 
relators (OjOk), requires the construction of the same 
number of transfer matrices (to cover the special case 
of the correlator where j = k), but these matrices must 
be multiplied and traced over up to separate times. 
Since some of these products of transfer matrices in these 
calculations will appear more than once, the actual com¬ 
putational cost can be reduced somewhat through use of 
a suitably “dynamically programmed” algorithm, where 
previously calculated products are saved and recycled [2] ■ 
Even in this case, however, far more than two solitary 
products would be required. Furthermore, as the or¬ 
der of the desired moment increases, the advantage 
of the moment-generating function method becomes in¬ 
creasingly pronounced, as the numerical derivative will 
require only approximately n expectation values, instead 
of L". 

Simply put, the fact that the exponential nature of the 
moment generating operator turns long Kronecker sums 
into simple Kronecker products makes it ideally suited 
for use with a matrix product state. In all cases, only 
a small number of expectation values must be computed 
in order to allow the calculation of a numerical deriva¬ 
tive, with each expectation value containing the operator 
at every site j. Moreover, application of these local 
operators does not increase the bond dimension of the 
state. 

Such moments can in principle also be computed by 
means of an MPO [13]. As a straightforward demonstra¬ 
tion of this, consider for example the MPO given by [52] 



h 0 


coupled with the boundary conditions 


(34) 


(<^l| = (0 1), (35) 


and 



(36) 


W={cI>l\1[CM 

i=i 

= Eo.- 

i=i 


so that 


(37) 

In this naive implementation, each application of W to 
\ip) increases the bond dimension of the state by a fac¬ 
tor of two, and thus, to calculate the moment in this 
way requires a bond dimension exponential in n. This in¬ 
crease can be overcome, if necessary, by means of a stan¬ 
dard truncation approximation, as done in the TEBD 
algorithm. Alternatively, a more sophisticated MPO can 
be constructed [SS] which represents the operator Af", 
but with a bond dimension of just n -I- 1, resulting in a 
procedure which still scales linearly. 


2. Sums of many-body operators 

We now examine the case of calculating the higher- 
order moments of a more general set of operators M = 
'^^OjOj+i...Oj+k- In other words, we consider opera¬ 
tors which are a sum of terms acting on at most k sites 
at a time. So long as k is finite, it remains possible to 
evaluate the moment generating functions with a single 
expectation value. This can be done by appealing to the 
same iTEBD technique [lU [54] widely used to simulate 
time evolution and imaginary time evolution in tensor 
network states. 

To begin, we partition the terms of the operator into 
classes of mutually commuting operators. This will re¬ 
quire at most k classes. For example, consider as an 
operator the two-body transverse Ising Hamiltonian 


H = ( 38 ) 

3 

For simplicity, let us write H in a manner which makes 
it explicitly a sum of two-body operators 

u = + + ( 39 ) 

3 

Then, the terms can be partitioned between the even 
and odd pairs of sites, and H can be written as 


To evaluate the moments, one defines the total MPO W 
to be 


77 — Heven a odd-! 


(40) 


A, 


• * • 


+ I {a] + a|+0 (41) 

jeven 

and 

HoM = - E + f (^1 + ^l+i) ■ (42) 

jodd 

The moment generating function therefore has the 
form F{a) = ^ admits a Suzuki-Trotter 

approximation |55j . To second order in a, this has the 
form 


„aH^. 


„+aH„ 




„aH, 


dd p 2 




(43) 


Higher order versions of the approximation have also 
been well-documented and can be easily substituted 
where greater precision is required [56]. 

Because and Ho^fi were explicitly constructed to 

be sums of mutually commuting terms, each exponential 
in the right hand side of Eq. (431 is now in the same 


Kronecker sum form as we had in the case of on-site op¬ 
erators, and each can therefore be equivalently expressed 
as a single tensor product of operations acting on the en¬ 
tire state at once, in the manner of Eq. (33). This is 
graphically depicted in Fig. Hence, by applying the 
three exponentials from Eq. ( |4^ in sequence, we can 
easily calculate the expectation value that represents the 
moment-generating function. 

From this point, evaluating the expectation value is 
no different than the case of on-site operators. The bulk 
of the numerical costs are therefore essentially the same 
for both on-site and many-body operator, with only one 
difference: applying these layers of exponential operators 
will increase the bond dimension of the system, and the 
size of these bonds may need to be “truncated” by some 
approximation scheme to keep the system numerically 
tractable. This, however, is a common practice in the 
field of MPS algorithms, and is easily done by means of 
a Schmidt decomposition (see for example HH). 

While our example considered an operator with fc = 2, 
that is, two-body interactions which could be partitioned 
into two internally commuting classes, the technique eas¬ 
ily generalizes to larger interactions. The Suzuki-Trotter 
approximations, for example, can be iteratively applied 
to an operator e^+^+c approximating 

in terms of and and then approximating . 

The ability to calculate moments and cumulants for 
an operator with many-body terms has a particularly 
useful application in the world of numerical state esti¬ 
mation. A common goal of tensor network algorithms 
is the calculation of an approximate numerical ground 
state, from a given Hamiltonian H. These algorithms 
are typically iterative in nature, gradually refining the 
approximation as the energy E tends towards Eq. It is 


a) 





FIG. 3. (a) Graphical representation of the moment¬ 
generating function F{a) = for an operator M = 

Since each term in M acts at only one site, the 
moment-generating operator possesses the same structure, 
even though the moments M” are fundamentally non-local, 
(b) The moment-generating function for an operator which 
is the sum of two-body terms and which possesses the form 
H = Hodd + He„en, such as the transverse Ising Hamiltonian 
defined in Eq. (38l. The operator is approximated by the 
second-order Suzuki-Trotter formula in Eq. (431, which pro¬ 
duces three “layers” of operations. Each layer is a sum of 
two-body terms. 


therefore often important to have a means of actively 
checking this convergence during the course of the algo¬ 
rithm. The variance (second cumulant) of the Hamilto¬ 
nian, (AH^) = (iJ^) — (H)'^ is well-suited to this task |2|. 
For e = there will be an exact eigenvalue Eex 

within e of the approximate energy E. In other words 


\E-Eex\<e. (44) 

This quantity can therefore be used as an upper bound 
on the convergence of the system, and gives a sufficient 
criterion for halting. Although epsilon is not guaranteed 
to be small as soon as the energy has converged, if we 
iterate our algorithms until epsilon is very small, we can 
be assured that the approximate ground state energy is 
very close to the true value Eq. 

In |2| , a dynamically programmed algorithm was given 
for computing (Ai/^) in the context of a finite matrix 
product state. The method presented here performs 
at least as efficiently in that case, with the advantage 
that it can also be applied to infinite (or indeed, higher¬ 
dimensional) systems. In |43j . the same quantity was pre¬ 
sented and evaluated by means of an MPO. As discussed 
above, the MPO technique can in principle be used as 
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an alternative method to compute other moments and 
cumulants as well, at the cost of allowing the bond di¬ 
mension to increase. 

For the case of the energy cumulant, our method is also 
particularly well-suited for use with the TEBD/iTEBD 
algorithms. We have previously remarked that for 
a Hamiltonian Ff, the calculation of the associated 
moment-generating operator is essentially identical 
to an imaginary time-evolution operator with a time step 
of 6 t = a. Each iteration of this algorithm therefore 
amounts to calculating from which the moment¬ 

generating function F{a) = (e“^) can easily be com¬ 
puted. If one also computes F{—a) = the error 

bound can therefore be computed very efficiently up to 
order in accordance with Eqs. (261 and (27). 

Performing this check at regular intervals throughout 
the evolution offers a halting condition to certify conver¬ 
gence of the energy. In some respects, this convergence 
criterion is superior to the typical methods, which often 
signal a halt when SE, the change in the approximate 
energy between two successive iteration steps, drops be¬ 
low some minimum value. Such a method can occasion¬ 
ally give a false sense of convergence when the algorithm 
“stalls out” and begins evolving only very slowly, despite 
remaining some distance from the ground state. The 
variance of the energy provides information not about 
the convergence of the algorithm but of the energy itself, 
by identifying when the system is very close to an exact 
eigenstate. However, in some cases care must be taken 
that the nearby eigenstate is in fact the ground state, 
and not some excitation. 


C. Evaluating generating functions on an infinite 
matrix product state 

At first glance, it may seem that moment-generating 
techniques discussed above cannot be applied to these 
systems, since the value of a quantity like M = Oj 
is clearly diverging, and only related limits like 


In other words, the desired quantity is simply the 
largest eigenvalue of the transfer matrix associated with 
the state and operator in question. Hence by defining 


h^{a) =logF^{a), (48) 

we find that we have access to the per-site limits of the 
cumulants even in the infinite case. In the same manner 
as with finite systems, they are given by the derivatives 
of with respect to a 


l''p\a) = lim —Kn + 0{a). (49) 

We note briefly some practical considerations that are 
important when evaluating Ip^ for a real matrix product 
state. First, algorithms for generating the states which 
are based on the iTEBD principle are likely to require the 
use of a two-site unit cell, even if the state is expected to 
possess only one-site translation invariance, as a result 
of the two-body nature of most parent Hamiltonians. In 
this case, of course, a two-site transfer matrix is required 
{£ = 2), and we must take a square root of its largest 
eigenvalue in order to recover the correct per-site limit. 

Additionally, we observe that for the second cumulant 
in particular, it can be particularly desirable to calculate 
using the characteristic function G = (e“^) instead of 
F. This is because one can then take advantage of the 
fact that 


lim = h^ia) + h^i-a) + 0{a'^) 

L—¥00 Lj 

= log Goo (a) + log Goo (-a) 

Then, combining the two log terms and observing that 
Goo (—a) = Gao{a)*, we have 


lim -K 2 = log (|Goo(a)P) + G(a^). (50) 

L—>-oo L 


- oo 

(M)= lim -^(O,) (45) 

3 

are well-defined. In this situation, however, while the 
moment-generating function F defined above may di¬ 
verge, one can still define and calculate the related quan¬ 
tity 


Foo = lim (6“*^)!/^. 

L—^oo 


(46) 


As discussed in section HI, quantities of this form can 
in fact be computed quite naturally. Using equation 21 
(see Fig. [^, clearly we have 


In other words, we can calculate the per-site second 
cumulant up to second order in a by evaluating Goo (a) 
only once, and without directly computing any numerical 
derivative. 


D. Higher-Dimensional States 

Although in this work we shall be applying these 
techniques to one-dimensional systems, there is noth¬ 
ing about the procedures above that cannot be immedi¬ 
ately generalized to finite-sized higher dimensional states. 
Consider for example a total magnetization-type opera¬ 
tor on a 2 dimensional, L x L square lattice, given by 


F^ 


lim = 

L—^oo 


X 

X 


max 


max 


(47) 


M = (51) 

3 k 
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whereOjfe represents a specific operator acting locally on 
site (j, k) of the lattice. Of course by representing both 
j and k by some composite index J (now running from 
1 to L^) we can immediately see that M is no different 
than the magnetization-like operators we considered in 
the one-dimensional case 


M = Y,Oj, (52) 

,7 


and hence that our earlier analysis goes through: the 
moment generating operator e“^ can still be written as 


eaM 



(53) 


This object is still an operator that acts only locally and 
whose expectation value can be evaluated all at once 
(graphically depicted in Fig. [^. 

If we consider instead an operator H which contains 
many-body terms (but for whom each term acts non- 
trivially only on a finite number of sites), one can play 
the same tricks as in one dimension, first partition¬ 
ing the terms into mutually commuting sets of terms 
Hi,H 2 + then expressing the moment-generating 

operator as 


and finally applying some form of Suzuki-Trotter approx¬ 
imation as described above to express the operator as a 
product of exponentials, each of which can be applied to 
the state all at once. For nearest-neighbor interaction 
terms on a square lattice, the procedure is essentially 
identical to the one dimensional case, except that one 
must use four classes instead of two: two for interactions 
in the horizontal direction, and two for the vertical. For 
operators with more complicated terms (such as a sum 
of “plaquettes,”) the number of partitions may be larger, 
but in general the costs do not grow rapidly despite the 
increase in system dimension. 

Hence, in either case, the associated moment¬ 
generating operators can be disentangled into tensor 
products or sequences of tensor products, even in higher 
dimensions. Simply put, the essential “power” of the 
moment-generating function method is the fact that the 
moment-generating operator of a local operator M 
is itself a local operator, and this fact is unchanged re¬ 
gardless of the dimensionality of the system. 

Once the moment-generating operator has been ex¬ 
pressed as a tensor product and applied to the state, it 
still remains to numerically contract the tensor network. 
It is at this stage where things become more difficult than 
in the one-dimensional case, since computing the expec¬ 
tation value of any operator on a higher-dimensional ten¬ 
sor network state can be quite hard. Exact calculation 



FIG. 4. The moment-generating operator for an 
operator of the form M = applied to a two- 

dimensional state on a square lattice. As in the one¬ 
dimensional case, the locality of each term in M ensures 
the locality of the terms in and hence, the moment¬ 

generating operator can still be evaluated all at once, by ap¬ 
plying the appropriate onsite operator at each lattice site. 


has been shown to be exponentially costly in L (in partic¬ 
ular, it is a ^P-hard problem m)- Nevertheless, a wide 
variety of numerical techniques have been developed to 
approximate these contractions efficiently with minimal 
errors, the details of which are outside the scope of this 
paper (See for example refs. [5HH55] ). We do caution, 
however, that in our experience, the higher the order of 
the moment, the higher the sensitivity of the result to 
the errors introduced by approximate contraction. 


V. EXAMPLES 

A. Spin-1/2 Transverse Ising Model 

As has become almost customary, we begin by demon¬ 
strating our technique in the context of the widely- 
studied transverse Ising model; a chain of spin-1/2 par¬ 
ticles governed by the Hamiltonian 

L 

H = ( 55 ) 

j 

This model is a useful proving ground as it has been 
extensively studied and admits a well-known analytical 
solution ES] , as well as possessing a straightforward 
order parameter of the total magnetization 

in the x-direction. We can therefore test our techniques 
by using them to study the phase transition known to 
occur exactly at He = 1. To apply the Binder cumulant 
technique, we first use a numerical method to find the 
ground state by solving the generalized eigenvalue prob¬ 
lem [5]. We consider system lengths between L = 10 and 
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L = 45 in steps of five, using a bond dimension of x = 10 
(we verify that increasing the bond dimension does not 
change the results of the methods up to our working pre¬ 
cision). For each system length, ground states of the 
Hamiltonian in Eq. (55) are computed as we sweep over 
a range of values for the held coefficient B. Then for each 
value of L and B, we compute the Binder cumulant using 
the methods described above, by hrst computing ^2 and 
/i 4 by means of the moment-generating function. 


As shown in Fig.[^ the crossings of the Binder cumu- 
lants at various lengths are already clustered very close 
to the transition point, even though the lengths of the 
states are relatively short compared to the thermody¬ 
namic limit. But the location of the critical point can 
be computed to even greater accuracy by considering the 
pattern of successive crossings. These crossings show a 
clear trend towards a limiting value as the system sizes 
increase. This limiting value can be estimated by means 
of the BST Algorithm [55], which has been found to be 
a very powerful tool for estimating the infinite limit of a 
series of data based on finite size corrections which obey 
a power-law, even from a relatively small number of data 
points m- From this extrapolation, we estimate a crit¬ 
ical point of Be = 1.001(1). Here and elsewhere, the 
reported uncertainty in our extrapolation represent an 
estimate of the typical effect of uncertainties in the loca¬ 
tion of the crossing points, propagated through the BST 
Algorithm. Much more detail about the BST technique 
can be found in Appendix A. 


open 



FIG. 5. (color online) A Binder cumulant study of the trans¬ 
verse Ising model. The cumulants are computed for different 
system sizes across a range of values for the transverse field 
B (some intermediate system sizes have been suppressed for 
clarity of the figure). Crossing points are interpolated for 
successive pairs of curves, i.e. L = 10 and L = 15. These 
crossing values can then be seen to approach the known value 
of the critical field. Be = 1 (inset). The BST algorithm is 
used to extrapolate these values to the infinite limit, which 
gives Be = 1.001(1). 


The critical point of the Ising model can be probed 
directly through the higher-order cumulants of the order 
parameter, as well. Through the techniques above, these 
can be calculated from the finite systems at various sys¬ 
tem sizes. Alternatively, we can calculate a ground state 
for the infinite system through the iTEBD algorithm (in 
this case using a bond dimension of x = 20) and then cal¬ 
culating the second cumulant directly. Both procedures 
are showcased in Fig.[^ The behavior in the infinite case 
can be seen to agree with the limiting trend of the finite 
systems as the length is increased. The cumulant can be 
seen to become singular near the critical point, and can 
also be used to detect the transition. Using this method, 
we estimate B^. = 1.00(1) 

This technique can also be applied to the magnetiza¬ 
tion in the transverse direction, In. this 

case, it is the derivative of the cumulant which becomes 
singular in the thermodynamic limit to signify the crit¬ 
ical point. Again, the results from the finite chains can 
be seen trending towards the inhnite limit (see Fig. [^. 



FIG. 6. (Color online) Per-site value of the second cumu¬ 
lant of the longitudinal magnetization, j-{AM^) = — 

{Mx)^), for the transverse Ising model. The cumulant is plot¬ 
ted for various finite system sizes, plotted against a range of 
applied fields. As the system length increases, the behavior 
tends towards the infinite limit (inset). In the limit, the cu¬ 
mulant diverges at the critical point. 

The critical exponent of the correlation length of the 
model can also be studied by means of the Binder cu¬ 
mulant. Once the critical point has been estimated, the 
curves of U 4 can then be plotted against L^/'^{B — Be) 
for various values of v. At the true critical exponent, the 
data should “collapse” to a single functional form inde¬ 
pendent of L, as seen in Fig. 

Finally, we can use this model to demonstrate the util¬ 
ity of the energy variance (AiJ^) in assessing numerical 
convergence, as described above. Starting from a random 
state, we apply the iTEBD algorithm with the transverse 
Ising Hamiltonian and evolve towards the ground state, 
over a range of field strengths B. As shown in Fig. 
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FIG. 7. (Color online) Second cumulant of the transverse 
magnetization, (AMf) = (Mf) — for the transverse 

Ising model (computed per site). The cumulant is plotted for 
various finite system sizes, plotted against a range of applied 
fields. As the system length increases, the behavior tends 
towards the infinite limit (inset). In the limit, the derivative 
of the cumulant shows a discontinuity at the critical point. 
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FIG. 8. (Color online) The Binder cumulants for the trans¬ 
verse Ising model, plotted for a variety of system sizes as a 
function of — Be) for the known values u = 1 and 

Sc = 1. As expected, for these values the curves are seen to 
collapse to a functional form essentially independent of the 
length scale. This property can be used to estimate the val¬ 
ues of the critical point and the critical exponent by treating 
them as fit parameters and optimizing the collapse. 


B. Spin-1 Transverse Ising Model 

We consider next the spin-1 generalization of the Ising 
model, with Hamiltonian 


h = -Y,s-s, 


i+1 


BS^. 


(56) 


Here, we have simply replaced the spin-1/2 Pauli ma¬ 
trices from Eq. (55) with their spin-1 counterparts. This 


model is of interest because unlike the spin- 1/2 case, 
it has no exact analytic solution. Nevertheless, in the 
thermodynamic limit the magnetization is qualitatively 
similar to the spin-1/2 case. Notably, it still displays a 
quantum phase transition at a critical value of the trans¬ 
verse field, which has been studied by various numeri¬ 
cal [5H1 US] . The accepted value for this critical field is 
Be = 1.326 [55]. 

We study this transition point with the same tech¬ 
niques as before, calculating the Binder cumulants from 
the second and fourth moments of the a;-magnetization 
for various system sizes (Fig. 10). The ground states are 


calculated using the same finite MPS technique and with 
a bond dimension of y = 20. The successive crossings 
are then compared and a limiting value extrapolated us¬ 
ing BTS. In this case, we compute an estimate of the 
transition at Be = 1.327(1). 

Once again, it is also constructive to consider the sec¬ 
ond cumulant on its own. Numerical calculations of the 
magnetization for this model invariably show some finite 
size effects around the transition, producing a finite “tail” 
near the transition point, which makes an exact determi¬ 
nation difficult using the order parameter alone. But 
the transition appears much more sharply as a singular¬ 
ity when we consider the second cumulant, as in Fig. |11| 
(higher cumulants such as K 4 can also be used for this 
purpose). From this quantity, we obtain Be = 1.324(2), 
an estimation which is to within less than 0 . 2 %. 

As before, can also study the critical exponent ly, 
known for this model to be the same as the spin- 1/2 case, 
1 / = 1 . As a proof of principle, the “data collapse” for 


the known values of v and Be are shown in Fig. 12 As 
expected, for these values the curves are seen to collapse 
to a functional form essentially independent of the length 
scale. This property can be used to estimate the values 
of the critical point and the critical exponent by treating 
them as fit parameters and optimizing the collapse. 


the results are initially somewhat noisy when compared 
to the analytically known E vs. B curve, which is re¬ 
flected by the large error bars computed from (Ai7^). 
However, as the algorithm continues, these error bars 
shrink and eventually become essentially zero, signaling 
the complete convergence of the energies. 


C. Spin-1 Ising Model in Crystal Field 

For another application, we consider also a variation on 
the spin-1 Ising model, where the usual transverse field 
has been replaced by a quadratic, crystal held term, to 
give the following Hamiltonian 
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FIG. 9. (Color online) The energy of the spin-1/2 transverse Ising model is calculated using approximate ground states 
generated by the iTEBD algorithm (x = 20). The numerical data (points) are plotted alongside the exact solution (line). Error 
bars are calculated from e = (Aff^), with (AJT^) the second cumulant of the energies, (a) After 10 steps, the energies are 

still noisy and the error bars are quite large, (b) After 20 steps, the error bars have clearly decreased, and are largest for the 
points with the largest discrepancies from the exact solution, (c) By 100 steps, the error bars are within the size of the data 
points, and the approximate energies are very close to the known analytical result. 



FIG. 10. (Color online) A Binder cumulant study of the 
spin-1 transverse Ising model. As above, the cumulants are 
computed for different system sizes across a range of values for 
the transverse held B (some intermediate system sizes have 
been suppressed for clarity of the hgure). Crossing points are 
interpolated for successive pairs of curves, i.e. L = 10 and 
L = 15, and the BST algorithm is used to extrapolate these 
values to the inhnite limit, which gives Be = 1.327(1). 
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FIG. 11. Per-site value of the second cumulant of the lon¬ 
gitudinal magnetization, — {Mx}^), com¬ 

puted for the Spin-1 Ising model. The cumulant is calculated 
for an inhnite system directly. 


L 

H = -J2S]S^^,+B{S]r. (57) 

3 

This variation of the spin-1 model admits a mapping to 
the spin-1/2 case (c.f. [70]), from which the critical point 
of i?c = 2 can be analytically obtained. To compare 
our method, we perform the same numerical calculations 
as above: first generating ground states at various finite 
lengths using MPS methods with a bond dimension of 
X = 10, and then computing the Binder cumulants. From 


the Binder curves (see Fig |13[ ) we once again perform the 
BST extrapolation of the successive crossings to arrive at 
an estimate of Be = 1.999(1). 


Direct examination of the second cumulant in the in¬ 
finite system is also still a viable method for estimating 
the transition point. In this case, the location of the 
maximum gives Be = 1.996(1) (see Fig. 141. As before, 
it also remains possible to study the critical exponent v 
and critical held value simultaneously by seeking to col¬ 
lapse the data to it s u niversal behavior as a function of 
L^I\B-Be) (Fig. [15 1. 
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FIG. 12. (Color online) The Binder cumulants for the spin-1 
transverse Ising model, plotted for a variety of system sizes 
as a function of — Be) for the known values zz = 1 and 

Be = 1.326. The length-independence of the curves allows 
this technique to be used as a means to estimate both v and 
Be 


D. Spin-1/2 Ising model on a 2D lattice 

Finally, as a proof-of-principle, we briefly demonstrate 
the application of these methods to a two-dimensional 
system: the spin-1/2 Ising model on a 2D square lattice. 
This system is described by the Hamiltonian 

^ = - E + BSlu, (58) 

i.fc 

where the subscripts are understood to terminate at 
the boundary of the system. As discussed above, the 
study of two dimensional systems with tensor network 
states is considerably more involved than the study of 
one-dimensional systems. More elaborate methods must 
be undertaken to numerically approximate the ground 
states, and elaborate approximation schemes must be 
used in order to calculate expectation values, which are 
otherwise prohibitively costly in time and memory. A 
detailed and high-precision study of the critical point of 
this model is therefore beyond the scope of this paper 
(see instead [H]). Nevertheless, we include the following 
rough estimation in order to demonstrate how easily the 
moment-generating function method can be generalized 
to higher dimensional states, as well as to underscore the 
utility of Binder cumulant techniques even for data cal¬ 
culated relatively cheaply. 

To this end, we generate approximate ground states 
for the model, using a simple method of local updates (a 
2-D generalization of TEBD) and the smallest nontriv¬ 
ial bond dimension, % = 2. To check the behavior, we 
consider the order parameter 



FIG. 13. (Color online) A Binder cumulant study of the 
spin-1 transverse Ising model with crystal field. As above, 
the cumulants are computed for different system sizes across 
a range of values for the transverse field B (some intermediate 
system sizes have been suppressed for clarity of the figure). 
Crossing points are interpolated for successive pairs of curves, 
i.e. L = 10 and L = 15, and the BST algorithm is used 
to extrapolate these values to the infinite limit, which gives 
Be = 1.999(1). 



B 

FIG. 14. Per-site value of the second cumulant of the longitu¬ 
dinal magnetization, 4(AM|) = — computed 

for the spin-1 Ising model with a transverse crystal field. The 
cumulant is calculated for an infinite system. 


3,k 




(59) 


Then, as in the one-dimensional case, we compute the 
Binder cumulant of the order parameter across a range 
of applied fields, for systems of size LxLuptoL=12, 
and observe the crossings (Fig. 161. The largest crossing 
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FIG. 15. (Color online) The Binder cumulants for the spin-1 
Ising model with crystal held, plotted for a variety of system 
sizes as a function of {B — Be) for the known values v = 1 
and Be = 2. As expected, for these values the curves are seen 
to collapse to a functional form essentially independent of the 
length scale. 


we are able to compute, between L = 10 and L = 12, 
occurs aX B = 3.11(1), which is already reasonably ac¬ 
curate compared to the accepted value of Be = 3.044, 
as calculated by quantum Monte Carlo [72]. A BST Ex¬ 
trapolation of the data gives Be = 3.3(2), a crude esti¬ 
mation but with relatively large error bars, owing largely 
to the fact that only three crossing values have been 
used in the extrapolation. We note also that, unlike the 
case of one-dimensional systems, the Binder crossings for 
two-dimensional systems are not necessarily converging 
monotonically and hence may not necessarily admit an 
easy extrapolation m- Instead, greater precision could 
likely be obtained through the use of more sophisticated 
two-dimensional methods (or additional computational 
resources) to study the crossings for slightly larger sys¬ 
tem sizes. 


VI. SUMMARY 

In this paper we have presented a method for efficiently 
calculating the higher order moments and cumulants of 
general operators on systems represented by tensor net¬ 
work states. For finite systems, this capability has a va¬ 
riety of applications in the search for phase transitions 
in quantum systems. Chief among these is the calcula¬ 
tion of the celebrated “Binder cumulant,” which provides 
a powerful tool for not only detecting phase transitions, 
but determining their location to a high degree of accu¬ 
racy using only relatively small finite systems to probe 
the infinite limit. The finite size scaling of the Binder cu¬ 
mulant also provides an estimate the critical exponent of 
the correlation length. Although the second cumulants 


FIG. 16. (Golor online) A rough Binder cumulant study of 
the spin-1/2 transverse Ising model on a square lattice, us¬ 
ing a local-update numerical algorithm with bond dimension 
X = 2. As in the one-dimensional case, the cumulants are 
computed for different system sizes across a range of values 
for the transverse field B. The largest crossing point, between 
L = 10 and L — 12, is at R = 3.11(1). Extrapolating to the 
infinite limit gives Be = 3.3(2), though this cannot be done 
with high reliability on such a limited dataset (see text). 


of Hamiltonians have been considered in the context of 
matrix product states, to our knowledge, critical point 
detection techniques based on the Binder cumulant (or 
cumulants in general) have not generally been put to use 
in studies based on tensor networks, despite being widely 
applied to classical systems and quantum Monte Carlo 
studies. It is our hope that the methods presented in 
this paper will allow them to be embraced by the tensor 
network community as well. 

In the case of infinite systems, we also present a 
method to calculate the per-site limits of the cumulants 
efficiently as well. The higher cumulants of an order pa¬ 
rameter often show sharp behavior at the critical points, 
which in many cases allows for easier detect than the 
changes in the order parameter itself. In particular, we 
show how singularities in the second cumulant can pro¬ 
duce a relatively precise (computationally cheap) estima¬ 
tion of the location of the transition. All the techniques 
(finite and infinite) are demonstrated in the context of 
the transverse Ising model, the spin-1 transverse Ising 
model, the Ising model in a crystal held, and could eas¬ 
ily be applied to other models. We also demonstrate a 
useful application of the second cumulant of the energy. 
This quantity, which we calculate for both hnite and in- 
hnite systems, provides a useful sufficient condition to 
determine when a numerical ground-state estimation al¬ 
gorithm has converged. As we demonstrate in the con¬ 
text of the Ising model, it can identify convergence up to 
a very high level of precision. 

Finally, we present a proof-of-principle demonstration 
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of the methods as applied to the transverse Ising model 
on a square lattice. Our result demonstrates that the 
methods for computing moments and cumulants are eas¬ 
ily generalized to states in two dimensions or higher. 
Precise calculation of the moments and cumulants of a 
such a system may be more difhcult, since state prepa¬ 
ration and the process of computing expectation values 
are themselves much more complicated in higher dimen¬ 
sions. However, the central idea of our method on its 
own remains just as straightforward as in one dimension. 

During the preparation of this work, it was brought to 
our attention that a technique based on matrix product 
operators (MPOs) has also been suggested as an alter¬ 
native method for evaluating cumulants and moments, 
particularly in the context of the second cumulant of the 
energy and its use as a convergence check. We believe 


both methods have complimentary strengths and weak¬ 
nesses, depending on the context in which they are ap¬ 
plied. 
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Appendix A: Extrapolation with the BST Algorithm 


We now briefly overview the Bulirsch-Stoer extrapola¬ 
tion scheme, commonly referred to as the “BST” Algo¬ 
rithm (the meaning of the “T” in this acronym has evi¬ 
dently been lost to time). This method was introduced in 
[66j in the context of differential equations, but has been 
widely adopted as an extrapolation scheme whenever one 
seeks to project a sequence of data with unknown func¬ 
tional form to its infinite limit. In particular it has be¬ 


come a useful tool for finite-size scaling techniques, and 
was studied extensively in the context of lattice models 

in [57] . 

The BST Algorithm assumes that we are attempting 
to extrapolate to a limiting value for an infinite system, 
which is subject to power-law corrections when approx¬ 
imated by a finite-size. For example suppose we have 
a sequence of critical field values which are approxi- 
mants to the true value of the critical field in the in¬ 
finite system: {B{Li), B{L 2 )...B{Ln)}, which approach 
Boo = B{L ^ oo). The BST Algorithm applies when, for 
each estimate B{L), B^o — B{L) = P{L) for some fixed 
(but unknown) polynomial P. This pattern of power-law 
corrections has generally been found to be true in the 
case of Binder Cumulants [TIES]. 

The technique works by taking the initial sequence 
02 °^ using it to construct a new se¬ 
quence, {a^\ a^\ whose convergence towards 

the infinite limit has been accelerated, so that is 

in fact a better estimate than . For clarity, note that 
we are using parenthetical superscripts to label the se¬ 
quence, and subscripts to enumerate the terms within a 
sequence. 

The terms in this new sequence are defined as follows 


„a+i) _ 


= a 
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Note that, since the denominator of equation ]Al| makes 
reference to the sequence it is necessary to define 

the sequence to handle the initial step of the 

algorithm in which j = 0. To that end, one simply takes 
= 0 for all k. 

This procedure can then be repeated, taking the se¬ 
quences and as the inputs to generate and 
so on. This iteration can be done at most N — 1 times, 
at which point the resulting sequence contains 

only one term. This term is the BST algorithm’s best 
estimation of the infinite limit of the original sequence. 

The parameter “w”, which appears as the exponent on 
the length scales, is a free parameter in the algorithm. 
The value of to which gives the best convergence will de¬ 
pend on the form of the power law corrections in the 
original sequence, which is generally unknown. Hence, in 
practice, a range of parameters must be considered, se¬ 
lecting the one which best optimizes the convergence. To 
this end, note that a rough estimate of the “precision” of 
the sequences can be made by 


(A2) 

This value should be decreasing with each iteration if 
the procedure is correctly accelerating the convergence 
of each new sequence. The final value of this estima¬ 
tor, A final = gives a convenient way to fix the 
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free parameter w: we simply repeat the algorithm while 
varying w, and choose the one which minimizes A final ■ 
In practice, it has often been observed ITS] that the de¬ 
pendence of the estimations on w is rather weak, with 
large ranges of values giving comparable results. In other 
words, it is often more important simply to avoid a “bad” 
value of uj than to try to find its absolute “best” value. In 
our work, we have used as our procedure a sweep over the 
range uj G (0,2], testing with steps of size 0.1. We also 
require that our extrapolation be ’’stable” under small 
variations in uj. 

We note that the value of A final cannot be used as a 
complete measure of the error in a final estimation. It is 
a measure of the internal consistency and the precision of 
the acceleration in the BST algorithm, but cannot con¬ 
tain any information about whether the algorithm has 
captured the ’’true” functional form of the finite-size cor¬ 
rections. Additionally, it does not reflect the propaga¬ 
tion of errors on the data which are being extrapolated. 


A very small and stable value of A final indicates that 
the algorithm is extrapolating the data to the best of it’s 
capability given its assumptions and the finite number 
of input points. It does not necessarily indicate that the 
result is extremely precise. 

In this paper, to estimate the error in a BST ex¬ 
trapolation, we start with the uncertainties of the input 
points, and essentially determine the propagated uncer¬ 
tainty empirically. In our case, the input points are the 
crossings of Binder cumulant curves. The crossings are 
computed by linearly extrapolating between data points, 
so we generously assume an uncertainty of one half the 
step size between points. The error is then estimated by 
considering a “worst-case scenario” in which the first few 
crossing points are perturbed downward by this amount, 
and the later points perturbed upwards. We run these 
perturbed points through the BST algorithm and observe 
the effect on the resulting extrapolation. The size of this 
effect is taken to be a rough upper bound on the total 
uncertainty. 


